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Canalization is a classic concept in Developmental Biology that is thought to be an important 
feature of evolving systems. In a Boolean network it is a form of network robustness in which a 
subset of the input signals control the behavior of a node regardless of the remaining input. It 
has been shown that Boolean networks can become canalized if they evolve through a frustrated 
competition between nodes. This was demonstrated for large networks in which each node had 
K = 3 inputs. Those networks evolve to a critical steady-state at the boarder of two phases of 
dynamical behavior. Moreover, the evolution of these networks was shown to be associated with 
the symmetry of the evolutionary dynamics. We extend these results to the more highly connected 
-R' > 3 cases and show that similar canalized critical steady states emerge with the same associated 
dynamical symmetry, but only if the evolutionary dynamics is biased toward homogeneous Boolean 
functions. 

PACS numbers: 89.75.Fb, 87.23. Kg, 05.65. +b, 89.75.Hc 



I. INTRODUCTION 

Boolean networks were originally proposed as models 
of genetic regulatory networks and are now widely used as 
models of self-regulatory behavior in biological, physical, 
social, and engineered systems They are designed 

to capture essential features of the complex dynamics of 
real networks by "coarse-graining" that assumes that dy- 
namical state of each node is Boolean, or simply on/off 
For example, in a model of a genetic regulatory sys- 
tem each node corresponds to a gene and its Boolean 
dynamical state refers to whether or not the gene is cur- 
rently being expressed. The regulatory interactions be- 
tween nodes are described by a directed graph in which 
the Boolean (output) state of each node is determined 
by a function of the states of the nodes connected to it 
with directed in-links. It has been shown that, despite 
their simplicity. Boolean networks capture many of the 
important features of the dynamics of real self-regulating 
networks, including biological genetic circuits 

Perhaps the most notable feature of Boolean networks 
is that they have two distinct phases of dynamical behav- 
ior. These two phases are called "frozen" and "chaotic" , 
and in random Boolean networks there is a continuous 
phase transition between them[lol-[l2|. The two phases 
can be distinguished by how a perturbation in the net- 
work spreads with time: in the frozen phase a perturba- 
tion decays with time, while in the chaotic phase a per- 
turbation grows with time [lol . [l3| . In networks in which 
the states of the nodes are updated synchronously, the 
two phases can also be distinguished by the distribution 
of network's attractor periods [3, 14]. When the updates 
are synchronous the system always settles onto a dynam- 
ical attractor of finite period. In the frozen phase the 
distribution of attractor periods is sharply peaked with 
a mean that is independent of the number of nodes N. 
In the chaotic phase the distribution of attractor periods 



is also sharply peaked, but with a mean that grows as 
exp(iV). In the "critical" state, at the boundary between 
the two phases, the distribution of attractor periods is 
broad, described by a power-law [l3 - [l7j . 

Many naturally occurring, as well as engineered, self- 
regulating network systems develop through some sort 
of evolutionary process. Motivated by this fact, a num- 
ber of models that evolve the structure and dynamics of 
Boolean networks have been studied p7l - l32j . These evo- 
lutionary Boolean network (EBN) models generally seek 
to determine the properties of networks that result from 
the evolutionary mechanism being considered. For exam- 
ple, some studies have explored evolutionary mechanisms 
that result in networks that have dynamics that are ro- 
bust again various types of perturbations, or that result 
in networks that are in a critical state. 

One example of an EBN is the model of competing 
Boolean nodes first introduced in Ref. |17j, and later 
studied in Refs. [2l[23,[2i|. In this model, the Boolean 
functions of the nodes evolve through a frustrated com- 
petition for limited resources between nodes that is a 
variant of the Minority game [33| . In the original pa- 
per on the model, it was shown that the network self- 
organizes to a nontrivial critical state with this evolu- 
tionary mechanism. Later it was discovered that this 
critical state is highly canalized [l^ . Canalization 34 [ is 



a type of network robustness, and is a classic idea in de- 
velopmental biology. Recently, experiments have demon- 
strated its existence in genetic regulatory networks (ssl - 
[37| . It occurs when certain expression states of only a 
subset of genes that regulate the expression of a particu- 
lar gene control the expression of that gene. Canalization 
is thought to be an important property of developmen- 
tal biological systems because it buffers their evolution, 
allowing greater underlying variation of the genome and 
its regulatory interactions before some deleterious varia- 
tion can be expressed phenotypically f38|. Critical net- 
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works are thought to important for gene regulation be- 
cause they can store and transfer more information than 
either frozen networks, which have a static output, or 
chaotic networks, which have a random output [2]. 

In the previous studies of the EBN with competing 
nodes, it was found that, for large networks in which 
each node has K = 3 in-links, the system evolves to a 
steady state that is both critical and highly canalized. 
The canalized nature of the evolved state was found by 
measuring the average frequency at which each of the 
256 possible Boolean functions of three variables occur. 
It was discovered that the functions organize into 14 dif- 
ferent classes in which all of the functions in each class 
occur with approximately the same frequency. Moreover, 
all functions in a class are equally canalizating and the 
classes whose functions are more canalizing occur with 
higher frequency. It was then found that the existence of 
these 14 classes is due to the symmetry properties of the 
evolutionary dynamics [39| . 

In this paper, we extend these results to more highly 
connected Boolean networks in which nodes have K > 3 
in-links. We show that, although an unbiased implemen- 
tation of the same evolutionary process does not lead to a 
critical state as it does for K = 3,hy biasing the process 
in a way that encourages the evolution of more homoge- 
neous functions the system can evolve to a critical steady 
state. We also show that the possible Boolean functions 
again organize into classes that depend on the symmetry 
of the evolutionary dynamics and that all functions in 
a class occur with the same average frequency. In this 
case however, the bias added to the evolutionary process 
causes a competition between the preference for canal- 
ization and a preference for homogeneity induced by the 
bias. This competition between canalization and homo- 
geneity is reflected in the frequency at which the Boolean 
functions occur. 



II. THE MODEL 

A. Random Boolean Networks 

A Random Boolean Network (RBN) consists of N 
nodes, i = 1, ....,N, each of which have a dynamical state 
with Boolean value ai = or 1. The Boolean state of 
each node is a function of the Boolean valued states of 
a set of Ki other randomly choosen nodes that regulate 
it. The regulatory interactions of the nodes are, thus, 
described by a directed graph. For synchronously updat- 
ing Boolean networks, which are the only ones considered 
here, the states Ci at time t -f 1 is a function of all states 
of its Ki regulatory nodes {ii, 12, i^. } at time t: 

a,{t + l) = M<j,,{t),a^,{t),...,aKAt)] (1) 

The function fi is a Boolean function of Ki inputs that 
determines the output of node i for all 2^* possible sets 
of input values. In this particular study, we consider 



RBNs in which Ki is fixed for all nodes i. No self-links, 
or multiple in-links from the same node are allowed in 
our models. 

In RBNs, each of the different functions fi is chosen 
randomly. We do this by choosing the 2^ outputs of each 
of the N functions either to be '0' with probability p and 
'1' with probability l—p, or to be '0' with probability 1—p 
and '1' with probability p. Which of the two cases is used 
is chosen with equal probability for each function fi , but 
remains fixed while assigning all of the individual outputs 
to that particular function. By choosing the Boolean 
functions in this way, a symmetry between 'O's and 'I's 
exists on average in the network. Note that the effective 
range of p is only from 0.5 to 1. 

Note that Hp = 1/2, the choice of functions is unbiased 
and each of the 2^^ functions are equally likely to be 
chosen. For p 1/2 the choice of functions is biased 
toward homogeneity. The homogeneity of a function fi 
is defined as the probability that it will output a "0", 
or the probability that it will receive a "1", whichever 
is larger, assuming that it receives random input. The 
average homogeneity of the network, P, is the average 
homogeneity over the functions of all nodes i. 

The system state E of the network at time t is given 
by the array of Boolean values of the states of each node: 

j:{t)^{ai{t),...,aN{t)} (2) 

Because the dynamics prescribed in Ec^. 1 are determin- 
istic, and the space of all possible network states is finite 
(of size 2^), all dynamical trajectories eventually become 
periodic. That is, after some possible transient behavior, 
each trajectory will repeat itself after some number of 
discrete time steps F to form a periodic cycle given by: 

I](i) = S(< + F) (3) 

The periodic trajectory over this cycle is referred to as 
the "attractor" of the dynamics, and the minimum F that 
satisfies this equation is the "period" of the attractor. 

As mentioned above, two distinct phases of dynamical 
behavior, "frozen" and "chaotic", exist for RBNs. For 
networks with uniform if > 2, networks are in the chaotic 
phase when p is near 0.5 and in the frozen or "fixed" 
phase when p is near or 1. There is a continuous phase 
transition between these phases at the so called "edge of 
chaos." The critical value Vr. where this transition occurs 
satisfies the expression [10l4l2 | 

l = 2Kp,il-p,) (4) 

We will refer to the networks that are critical because 
they satisfy Eq. 4 by construction as "non-evolutionary" 
networks, since such networks have no evolutionary dy- 
namics associated with them. 



B. Evolutionary Game of Competing Nodes 

The process for evolving the Boolean functions in the 
network is: 
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1. Start with an RBN constructed with a bias p, and 
choose a random initial state S(0). 

2. Update the state of the network using Eq. 1, and 
determine the attractor of the dynamics. (The at- 
tractor can be found using the algorithm discussed 
in the Appendix of Ref. |25l) 

3. For each update on the attractor, determine which 
Boolean value is the output state of the majority 
of the nodes, and give a point to each node that is 
part of the majority. 

4. Determine the node with the largest number of 
points; that node "loses". If two, or more, nodes 
are tied, pick one at random to be the loser. 

5. Replace the function of the losing node with a new 
randomly chosen Boolean function with bias p. 

6. Return to step 2. 

The essential features of the game are (1) frustration 
[sl,!!^!, since most nodes lose each time step, (2) negative 
reinforcement, since losing behavior is punished, and (3) 
extremal dynamics [4l|, since only the worst performing 
nodes Boolean function is changed. 

Note that previous studies of this evolutionary model 
have always replaced the function of the losing node with 
an unbiased, p — 0.5 random Boolean function. Thus, 
the evolutionary game we study here differs from previous 
work only at step 5. 

If the attractor period is longer than some limiting 
time, Fmax, then the score is kept only over that limited 
time. In our simulations, Fmax — was used. Each 
progression through the game is called an "epoch" . 

In this evolutionary process only the Boolean functions 
evolve; the directed network describing the regulatory in- 
teractions between nodes does not change. However, as 
we will see, canalization effectively removes interactions 
from the network. Thus, as the Boolean functions evolve 
and canalization increases the directed network of regu- 
latory interactions effectively changes j39j. Thus, effec- 
tively, both the structure of the network and the dynam- 
ics of the nodes simultaneously evolve, making the model 
effectively a co-evolving adaptive network model >42] . 

III. CANALIZATION AND SYMMETRY 

A. Canalization and Ising Hypercubes 

As mentioned above, for large networks with K = 3, 
it has been found [H ^ that the 256 = 2^^^ possi- 
ble Boolean functions of three variables organize into 14 
different classes in which the functions belonging to each 
class occur with the same frequency in the critical steady- 
state that results from the evolutionary game. Moreover, 
all the functions in a class are equally canalizing, and, 
in the steady state, functions in classes that are more 



canalizing occur with higher frequency. Thus, for K = 3 
the game causes networks to evolve to a critical steady 
state that is highly canalized. Canalization occurs in 
Boolean networks when the Boolean functions assigned 
to the nodes are canalizing. A Boolean function is canal- 
izing if its output is fully determined by a specific value 
of one, or more of its inputs, regardless of the value of 
the other inputs. The canalization of a function can be 
quantified by a set of numbers Pfc,fc = 0,1, — 1, 
which are defined as the fraction of the different possible 
sets of k input values that are canalizing [2]- 

Canalization can be further understood by mapping 
Boolean functions of K inputs onto configurations, or 
"colorings", of the ii'-dimensional Ising hypercube (39| . 
The Ising hypercube is a hypercube which has each vertex 
labeled, or colored, either '0' or '1' ("black" or "white"). 
In this representation, each of the 2^ possible sets of 
input values corresponds to coordinates on a given axis 
of the hypercube. The color of each vertex represents 
the output value of the function for the associated input 
values. The mapping of a Boolean function of K inputs 
to a configuration of the i^T-dimensional hypercube is one- 
to-one. 

This representation of Boolean functions as colorings of 
Ising hypercubes facilitates quick recognition of canaliz- 
ing functions. For a if-dimensional hypercube, the frac- 
tions of canalizing inputs of a Boolean functions are 
the fraction of its K — k dimensional hypersurfaces that 
are homogeneously colored. 

B. Symmetry of Evolutionary Dynamics 

This Ising hypercube representation is particularly 
helpful for understanding the role of symmetry in the evo- 
lutionary dynamics. The fourteen (14) different classes 
that were observed empirically in the original K = 3 
study were later recognized as being those Ising hyper- 
cube colorings that are related by cubic symmetry plus 
parity j39| . Parity in this case refers to simultaneously in- 
verting the Boolean values associated with each vertices. 
An underlying symmetry of the evolutionary dynamics 
was therefore reflected in this symmetry of the evolved 
steady-state. Furthermore, this symmetry preserves the 
canalization values Pfc of functions in each class, since 
neither cubic nor parity operations change the percent- 
age of homogeneously colored hypersurfaces. 

In mathematical terms the different classes correspond 
to the group orbits of the "Zyklenzeiger" group, which is 
the hyper-octahedral symmetry group 0„ (where n = K) 
combined with parity [39] . A group orbit is the set of con- 
figurations that map into each other through applications 
of a group's symmetry operations. The number of orbits 
can be calculated analytically using Polya's theorem 



where G is the symmetry group acting on the K- 
dimensional Ising hypercube, \G\ is the number of op- 
erators g E G, X is the set of hypercube colorings, is 
the set of colorings that are left invariant by g, and \X^\ 
is the size of set X^ . 

To apply this theorem, one must construct all the op- 
erators of the group, sum the number of functions left 
invariant by each operator, and divide by the total num- 
ber of symmetry operators. The hyper-octahedral group 
On has n!2" operations. Including parity operations dou- 
bles this number. 

A given symmetry operator g can be written as a per- 
mutation of the vertex numbers on a given hypercube. 
As a result, each operator g can be expressed in terms 
of its cycle structure x^2 ■■■^^in ^ where X^I^lo ~ '^^ ■ 
This notation indicates that g contains bi cycles of length 
1, 62 cycles of length 2, etc. The complete cycle repre- 



sentation of the hype-octahedral group for an arbi tra.; 
dimension K is given by a known recursion relation [4J] 
For K = the complete cycle relation is: 



5000 10000 

Epochs 

FIG. 1: (Color online) A graph of attractor period vs. epochs 
for a simulation of a = 4 network realization with an evo- 
lutionary bias of p = 0.65. Notice that neither short periods 
nor long periods dominate the behavior after the steady state 
is reached at ~ 10^ epochs. 
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IV. RESULTS 



where the coefficients of this polynomial indicate the 
number operators with a particular cycle structure. 

Without parity, the number of functions left invariant 
is equal to 2^"^, where iVc — SilLi the total num- 
ber of cycles in the operator. Parity must be treated 
separately; no functions are left invariant by the parity 
operator with any hyper-octahedral operator containing 
at least one cycle of length 1. Thus there are 2^^ func- 
tions left invariant for the operators which include parity, 
where Np = (1 — Q{bi)) X^t^i ^nd Q is the Heaviside 
step function. 

Thus, applying Polya's theorem to the = 3 case, we 
arrive at 

Pa = (l/96)((2)8 + 13(2)4 + 13(2)^ + 8(2)^ + 
8(2)2 _^ g^2)2 + Q{2f + 12(2)2 ^ 12(2)2) ^ ^j-^ 

This is precisely how many function classes were found 
empirically in the K = 3 case [1^ [s^. Similarly, the 
complete cycle representation for the if = 4 is 
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12X1 Xn '~\~ 51Xn 



12xj^X2 ~l~ 32x-|^X2 ~|~ 
4Sxlx2xl + 84xi + 96x1x1 + 48x1 



(8) 



Using this result, accounting for parity, and apply- 
ing Polya's theorem we calculate Pq = 222 for a 4- 
dimensional hypercube under rotation plus parity sym- 
metry. (Note that an erroneous value of Pq — 238 was re- 
ported in Rcf. 39] .) This is how many classes of functions 
should be observed in a critical steady state oi K — A net- 
works undergoing the evolutionary dynamics. Below we 
show that results from numerical simulations are consis- 
tent with this prediction. Moreover, the frequencies that 
these functions occur show a preference for canalization. 



A. Critical States of K=4 Networks with 
Competing Nodes 

We have performed exhaustive simulations of ensem- 
bles of networks with K — A playing the game of com- 
peting Boolean nodes. All of the simulation results re- 
ported in this subsection are for networks with N = 999 
nodes. Simulations were run for evolutionary processes 
with biases oi p — 0.5, p — 0.65, and p — 0.75. The 
frequency that each of the 2^ — 65536 different func- 
tions occurred in the evolved steady state was measured 
by simulating, for each p, an ensemble of 13,000 indepen- 
dent network realizations. Each realization was initial- 
ized with an independent random network with random 
links and different random functions biased with the as- 
sociated p value. The simulation of each realization was 
run for 10'* epochs to allow the network to reach a steady 
state. At the end the simulation of each realization, the 
functions of each node were recorded and then used to 
calculate the average frequency of each function for the 
ensemble of realizations. 

Figure 1 shows a graph of the attractor period vs. 
epoch for a simulation of a if = 4 network realization 
with evolutionary bias p — 0.65. For the first 6834 
epochs all the attractors found have periods longer than 
Tmax = 10"*. Then, attractors with shorter periods begin 
to appear. After about 10* epochs the network reaches 
a steady state with a broad distribution of attractor pe- 
riods. As the figure indicates neither chaotic behavior, 
exemplified by almost entirely large attractor periods, 
nor frozen behavior, exemplified by almost entirely short 
periods, dominate the behavior of the steady state. The 
steady state is instead at the "edge of chaos" and is a 
critical state. All network realizations for p = 0.65 and 
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FIG. 2: (Color online) Ensemble averaged frequency of K = 4 
Boolean functions in the evolved critical steady state. The 
black solid line corresponds to p = 0.65, and the red dotted 
line corresponds to p = 0.75. Functions are grouped together 
in classes based on hyper-octahedral plus parity symmetry 
and with groups then ordered according to canalization from 
high to low. The inset shows the same data for functions 
1-5000. 

p = 0.75 showed similar behavior once reaching a steady 
state. However, for p — 0.5, no attractor periods of less 
than 10,000 were observed for any network realization. 
That is, unlike in the K = 3 case, the system does not 
self-organize or evolve into a critical state when removing 
functions and exchanging them with unbiased functions. 
This indicates that the increased complexity of the more 
highly connected network disrupts the networks ability 
to self-organize with unbiased functions. 

These results indicate that a phase transition occurs in 
the evolutionary dynamics as p is increased. At values of 
p below the transition value the evolutionary dynamics 
do not produce a critical steady state, while at values of 
p above the transition, a critical steady state evolves. In 
fact, a second transition occurs at higher values of p. At 
values of p above this second transition, the evolved state 
is no longer critical but instead remains in the frozen 
state. Thus, critical steady states evolve only when the 
bias of the evolutionary process is within a range. 

Note that non-evolutionary K = 4 RBNs are con- 
structed with a critical bias of Pc ~ 0.85355, accord- 
ing to Eq. 4. This value of Pc produces networks with 
an average homogeneity of all Boolean function in the 
network P w 0.85358. However, the steady state value 
of P for networks evolved with a bias of p = 0.65 is 
« 0.71 Clearly, the critical state oi K = A evolved net- 
works is significantly different than the critical state of 
non-evolutionary K — A RBNs. 

Figure 2 shows the ensemble averaged frequency at 
which each of the 2^ — 65536 different K = A Boolean 
functions occur in the evolved steady state, for bias pa- 
rameters p — 0.65 and p — 0.75. The Boolean functions 
were ordered by first grouping them by their membership 
to a particular class under hyper-octahedral plus parity 
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FIG. 3: (Color online) Enlargement of a small portion of 
Fig. 2. The two classes of functions shown have the same ho- 
mogeneity but different canalization. Notice that for p — 0.65 
(black solid line), the class on the left with higher canalization 
occurs with a higher frequency than the class on the right with 
lower canalization. However, for p — 0.75 (red dotted line) 
both classes occur with approximately the same frequency. 
This indicates that the bias towards homogeneity is dominat- 
ing the behavior over the drive for canalization at the high 
value of bias for these two classes. 



symmetry. Then the classes were ordered in descending 
order by how canalizing the functions in the class are, as 
measured by the sum of their Pfe values. As expected, the 
graph shows that functions belonging to the same class 
occur with the same probability, at least to within the 
resolution allowed by statistical fluctuations. This con- 
firms the hypothesis that the underlying symmetry of the 
evolutionary dynamics is hyper-octahedral plus parity. 

Clearly, certain classes of functions occur with a higher 
probability than others. In general, functions on the 
left side of the graph (higher canalization) occur much 
more frequently than functions on the right (lower canal- 
ization). However, unlike in the K = 3 case, certain 
function classes with higher canalization occur less fre- 
quently than function classes with lower canalization. 
This occurs because, unlike in the previous K — 3 studies 
[itI [23I [23 . [26| , the evolutionary process here is biased 
toward homogeneity. In this case, the drive for canaliza- 
tion caused by the evolutionary dynamics is competing 
against the bias toward homogeneity for certain classes 
of functions. See Fig. 3. Nonetheless, the evolved critical 
steady state of these biased networks still shows a pref- 
erence for canalization and is strikingly different than a 
critical state of a non-evolutionary RBN where homo- 
geneity entirely dominates the relative frequency of func- 
tions. 

Analogous results presumably hold for even larger val- 
ues of K. However, because the number of Boolean func- 
tions for a given K goes as 2^ , accurately measuring 
the frequency that each function occurs at in the critical 
steady state becomes unfeasable for values of K greater 
than 4. 
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FIG. 4: (Color online) Critical state evolution order parame- 
ter as a function of evolution bias parameter p for networks 
of size N = 999 with K = 2, (black straight line), K ^ A (red 
dotted line), K = h (green dashed line), and K = % (blue 
dashed-dotted line). 



FIG. 5: (Color online) p when ^ = 0.5 as a function of 1/K 
for K = 4,5, 8. The straight (red) line is a linear fit to these 
three data points. 



B. Criticality as a function of p and K 

Our results for K — 4 indicate that canalized critical 
states evolve only when the bias p of the evolutionary 
game is within a range. At too low a value of p, the 
evolved steady state is a chaotic state and only long at- 
tractor periods are found. At too high a value of p, the 
evolved steady state is a frozen state and only short at- 
tractor periods are found. Only in an intermediate range 
of p does a critical steady state evolve. In order to quan- 
titatively find the approximate the range of p for which 
critical steady states evolve, we define an order parame- 
ter ^ as the percentage of steady state attractor periods 
that occur that are less than Tmax- Then, when ^ = 
the networks are assumed to be in the chaotic state, when 
^' = 1 the networks are assumed to be in the frozen state, 
and when < '5 < 1 the networks are assumed to be in 
a critical state. 

Figure 4 shows a graph of the order parameter 5" as 
a function of the evolution bias parameter p for network 
connectivities K = 3,4, 5, and 8. This data was produced 
using networks of size N = 999, using an equilibration 
time to reach the steady state of 10^ epochs, and then 
computing 4" over 3 x 10^ epochs. The ^ values were 
also then averaged over 140 network realizations for each 
value of p and K. (Sixteen realizations were used for 
K = 8). A period cutoff value of Tmax — lO"* was used 
in these simulations. 

For K ~ 3 the network already evolves to a critical 
state at p = 0.5, the smallest possible value of p, and 
stops evolving to a critical state at p sa 0.82. This is con- 
sistent with previous results [17| that unbiased K ~ 3 
networks evolve to critical states. This is not the case, 
however, for networks with K > 3. As shown in Fig. 4, 
the onset of evolution to criticality for these more highly 
connected networks occurs at a value p > 0.5. At least 



for K — 4, 5, and 8, and presumably for all finite values 
of K, evolution to a critical state occurs only for a finite 
range of p. For example, for K = 4, this range is from 
p « 0.60 to p oi 0.83. The width of this range appears to 
decrease, and both the minimum and median values of p 
appear to increase, as K increases. The findings from the 
results shown in Fig. 4 vary quantitatively, but remain 
qualitatively consistent if 1) the equilibration time is in- 
creased, 2) the value of Tmax is varied, or 3) the number 
of nodes N nodes is changed. 

Figure 5 shows the value of p when ^' = 0.5, which 
approximates the median value of p in the range of the 
evolution of a critical state, as a function of 1/K for 
K = 4, 5, and 8. Unfortunately, simulating networks 
with K ^ 8 is computationally unfeasable with our 
methods, and we are thus restricted to predict the asymp- 
totic behavior of these evolutionary random Boolean net- 
works using these relatively small values of K. The three 
points fall roughly on a straight line. If we extrapolate 
the linear fit of the data points, the value of p tends to- 
ward a value slightly larger than 1 in the limit of large K. 
However, this is physically unrealizable since p cannot be 
larger than 1. Therefore, we expect that as if — )■ oo, the 
width of the range of p for which criticality occurs goes 
to 0, while the median value of the range goes to 1. 

It is important to note that the range of evolution bias 
parameter p for which critical state evolves is, at all stud- 
ied values of K, less than the critical bias value Pc for non- 
evolutionary RBNs given by Eq. 1. Therefore, the critical 
steady state that results from evolutionary process is dif- 
ferent than the critical state of non-evolutionary RBNs. 
From previous studies of i^T = 3 networks, and from the 
results shown in Fig. 2 that were discussed above, the 
difference is that the evolved critical steady states are 
more canalized. 
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V. DISCUSSION 

In this paper we have extended previous work on an 
EBN model in which the nodes compete in a frustrated 
game that causes the Boolean functions of the nodes to 
evolve. The previous studies of this model found that 
the game causes the system to evolve to a critical steady 
state that is highly canalized. Canalized states and their 
evolutionary mechanisms are important in Developmen- 
tal Biology beause of the usefulness of the robustness 
against diliterous phenotypic expression that canaliza- 
tion in the genome provides. The previous studies also 
found that the evolutionary dynamics of the K = 3 
model has the symmetry of the 3-dimensional Zyklen- 
zeiger group, which is the combination of parity and the 
cubic symmetry group. 

The previous studies of this EBN, however, only con- 
sidered networks with K = 3 in-links per node. Here we 
extended the study the more highly connected networks 
with larger K. Real self-regulatory systems, both biolog- 
ical and engineered, typically have nodes with a range of 
K inputs Thus, it is important to understand how 
larger regulatory connectivity effects evolutionary mech- 
anisms. 

For networks with .ftT > 3, we find that the game as 
previously studied does not cause the network to evolve 
to a critical steady state. Instead, it will evolve to a 
chaotic steady state. This occurs because the unbiased 
game, which was studied previously, replaces the Boolean 
function of nodes that lose the game with randomly se- 
lected new functions that are choosen unbiasedly from 



the set of all possible Boolean functions. Apparently, for 
networks with K > 3, unlike what happens for networks 
with K = 3, a they are composed largely of nodes with 
random Boolean functions with an unbiased distribution, 
then the evolutionary game is not "strong enough" to in- 
duce a shift in the distribution of the nodes' Boolean 
functions sufficient to have critical state dynamics. 

However, we have shown that for networks with K > 3 
if the game replaces the Boolean function of the losing 
nodes with functions biased toward homogeneity, then a 
critical steady state can evolve. We studied the range 
of evolutionary bias that will cause critical state evolu- 
tion and found that it narrows and that its median in- 
creases with K. We have also shown that the critical 
steady states that evolve for K > 3 are highly canalized, 
although there is also a competing bias toward more ho- 
mogeneous Boolean functions. All functions in an orbit 
of the iiT-dimensional Zyklenzeiger group have both equal 
canalization and equal homogeneity and occur with equal 
frequency in the steady state. Thus, the symmetry of the 
evolutionary dynamics of the EBN with K regulatory 
links per node is that of the if -dimensional Zyklenzeiger 
group. 

This study illustrates the importance of symmetry in 
self-regulatory networks and of evolutionary processes. 
It would be interesting to analyze other self-regulatory 
network systems, both real and model systems, with the 
methods we have used. This would allow the importance 
of symmetry in evolutionary processes to be understood 
and become better appreciated. 
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